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Abstract 

We investigate the emergence of large, localized, pseudo-stable configurations (oscil- 
lons) from inflaton fragmentation at the end of inflation. We predict the number density 
of large oscillons, and the conditions necessary for their emergence in a class of inflation- 
ary models. Analytic estimates are provided for a 3 -|- 1 and 1 -|- 1-dimensional universe. 
We test our predictions with detailed numerical simulations in 1 -|- 1-dimensions. We 
see a zoo of oscillons emerging from the simulations, including the usual small ampli- 
tude "sech" oscillons as well as large "flat-topped" oscillons. The emergent oscillons 
account for ~ 80 per cent of the energy density of the inflaton. 
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1 Introduction 



Inflation [H El |3] is a wonderful mechanism for generating large scale density fluctuations 
that, under the influence of gravity, eventually result in the formation of structure in the late 
universe. However, inflation must eventually end to give rise to a hot, radiation dominated 
universe consistent with the success of Big Bang Nucleosynthesis (BBN) jHEj. Hence, we re- 
quire the inflaton to eventually decay into Standard Model fields, possibly via intermediaries 

Before the universe thermalizes, the inflaton often undergoes complex spatio-temoporal 
dynamics. In many cases it fragments on time scales t <^ H~^, leading to a turbulent, 
incoherent state of scalar waves (see for example [IDlIIIlIiailSlElIISlIISlinillSlIISlEniEIl 
[22] and for recent reviews and further references see [23l[2l])- In this paper, we discuss a class 
of real- valued, single field models, where the fragmentation can lead to copious formations of 
remarkably long-lived {t ^ H^^), localized, pseudo-solitonic configurations called oscillons. 

Oscillons have been known to exist in non-linear scalar field theories for some time. The 
earhest investigations can be found in ESI [2Z! (also see[2Hll2SlEniEIlE21E3l[MlESlEniEZl 
EHlEHlllOlllIlliailalllilllSllle]). They are localized in space and oscillatory in timeQ They 
are long lived, surviving for many thousands of oscillations. They are quasi-stable, losing 
energy very slowly through outgoing scalar radiation (for example, see [H], HH HHl HHl 150]). 
Unlike Q-balls j5Tl|52], these configurations arise in real fields and have no conserved charge, 
although see [HS] for an adiabatic invariant. Their quantum mechanical decay rate was 
recently investigated in t46||^ 

The presence of oscillons could have important consequences for the post inflationary 
universe. They could lead to an enhanced decay rate for the inflaton [51] or explosive 
production of particles in localized regions [Ml B6] , long after the homogenous inflaton has 
fragmented. Enhanced density perturbations due to oscillons provide seeds for formation of 
structure via accretion and mergers on much smaller scales than those responsible for current 
large scale structure. For large oscillons, there is a possibility of producing primodial black 
holes (for example, see [55]), which could in turn constrain the inflationary potential. In 
addition, it will be interesting to study gravitational wave production due to the formation, 
collapse and mergers of these lumps. A related investigation regarding gravitational waves 
from gravitationally collapsed inflaton clumps was carried out recently in [56j and from 
Q— balls in [53, EH] |59l EO]. In models where the field fragments efficiently into oscillons, the 
gravitational wave spectrum from preheating [HH [62l [631 EH [HHl EE] could be modified. None 
of these consequences of oscillons have been properly explored in a cosmological setting. 

An important first step towards exploring the cosmological consequences of oscillons 
is estimating their number density, fraction of energy density in oscillons, and individual 

* Hence they are not constrained by Derrick's theorem [4^ which only apphes to static configurations 
^They eventually decay either due to couplings to other field or their own outgoing radiation. The 
timescales depend on the model under consideration. 
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characteristics. In this paper, we undertake this task for a class of single field inflation 
models (or effectively single field models during the oscillatory phase). Near the minimum 
of the potential, the dependence on the inflaton is assumed to be V{ip) = rn^ip^ /2 — 
A<^^/4 + (7^m~^(y9^/6 . . .. We require (i) V'{(p) — m'^ip < for some range of the field and 
(ii) (X/g)"^ ^ 1. Condition (i) is generic for any potential that flattens out at large field 
values, for example the potential in axion monodromy inflation [UT] or models of hybrid 
inflation [BS]. Note that (i) does not imply the need for an inflection point, whereas (ii) 
helps in protecting the oscillons from a collapse instability, though is not strictly required 
[15] . Throughout the analysis, we assume that the self interactions of the inflaton dominates 
over the coupling to other fields during the time of interest. We discussed this model in 
detail in including an analysis of oscillon solutions in an expanding universe, and their 
stability. 

We now provide a brief synopsis of our approach and the organization of the rest of 
the paper. Details about our choice of potential and associated assumptions are discussed 
in section [2j In section |3] we start with zero point fluctuations in the inflaton during its 
oscillatory phase and follow their linear evolution analytically. Fluctuations in a limited band 
of wavenumbers get amplified via parametric resonance. Non-linear effects that give rise to 
large oscillons only become important if the initial fluctuations are amplified sufficiently fast, 
Hubble expansion being the competing effect. We use this criterion to identify the values of 
A, g and m for which we get copious production of oscillons. This linear analysis also allows 
us to determine an important scale, k^i, the scale that is the first to become nonlinear. 
We hypothesize that it is this scale that determines the co-moving number density of large 
oscillons: noscO'^ ~ {k^i/2TT)^. In section |4] we obtain fcni and hence estimate the number 
density in terms of the parameters of the inflaton Lagrangian (m, A and g). 

To test whether our estimate provides a good approximation for the number density, we 
need a lattice simulation of the field in an expanding universe. In this paper (section |5]) 
we discuss a detailed set of numerical simulations in 1 -|- 1-dimensions. We investigate the 
individual characteristics of oscillons, the fraction of energy density in oscillons, as well the 
number density of oscillons produced and compare them to our analytic estimates. We find 
good agreement between the two (better than a factor of 2) and the expected scaling as 
parameter values are changed. Results from numerical simulations in 3 + 1-dimensions will 
be presented in an upcoming publication [SHlj^ 

Our conclusions and future directions are presented in section |6| The details of our 
numerical set-up, additional discussion of initial conditions and a derivation of the Floquet 
exponent are deferred to appendices. Animations based on our numerical simulations can 
be found online at http://www.mit.edu/~mamin/oscillons.html, 

"''Given the range of scales that need to be resolved {H < k/a < m where H <^ m), the 3 + 1 dimensional 
simulations are quite challenging. Preliminary simulations in 3 + 1-dimensions (run in collaboration with 
Richard Easther and Hal Finkel) are consistent with the ansatz presented in this paper to within a factor of 
~ 2 ^6^. In 3 + 1-dimensions, oscillons also take up > 50% of the energy density of the inflaton, as is the 
case in 1 -|- 1-dimensions. 
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Before we begin our investigation, we briefly review the relevant literature on emergence 
of oscillons in the early universe. This list is by no means exhaustive. Formation of "axitons" 
in the axion field during the QCD phase transition was investigated in [70]. Emergence of 
such pseudo-solitons (not all of them are oscillons) has been investigated in certain super- 
symmetric, hybrid inflation models [Ml [TTj [72] . The approach in [51] is closest to the one 
taken in this paper. However, the potential considered in [HI] makes the oscillons susceptible 
to a collapse instability in 3 + 1-dimensions (see [IS])- Our method for calculating kj^i is some- 
what different, specifically revealing the impact of the expansion rate on the number density. 
Q-ball production in the early universe shares many similarities with oscillon production 
[571 EH EH] , except, unlike oscillons, Q balls are present in complex valued scalar field with 
an exactly conserved charge. Formation of oscillons from quasi-thermal initial conditions in 
a 1 + 1 dimensional de-Sitter universe was investigated in [SO]- A numerical investigation of 
oscillon-oscillon and oscillon-domain wall interactions in 2 -f- 1 dimensions was carried out in 
[38] . Very recently, results of a 3 + 1 dimensional numerical simulation of oscillons emerging 
from symmetry breaking phase transitions (with quasi-thermal initial conditions in deSitter 
space) was presented in [73]. They reported ~ 2% of the total energy density fraction in 
oscillons. To the best of our knowledge, an estimate for the number density in terms of the 
parameters of the model and its comparison with simulations was not provided in any of the 
above mentioned papers. 



2 Our model and associated assumptions 

We begin with the an action which includes the inflaton and gravity {h = c = 1) 



s = I dxW-g 



T. 



where Q is the determinant of the metric, R is the Ricci scalar and rripi is the reduced Planck 
mass. We assume that 

V(^) = lmV--,^' + £,y + ^^- (2) 

near the minimum of the potential (see figure [T]). The effective mass of the inflaton at the 
bottom of the potential is m whereas A and g are dimensionless parameters. We will assume 
that \X\,g ^ 1 and m/mp\ <^ 1. The fiducial values we have in mind are m/m^i ~ 5 x 10~^ 
and |A|,(7 ~ 10^^. The typical field values of interest are <^ m/y/X. For such values, we 
assume that we can ignore terms beyond ip^. Furthermore, we assume that far away from 
the minimum, the potential is consistent with results from the cosmic microwave background 
anisotropies [71]. We note that although we assume m/rripi ^ 5 x 10~^ (chosen to be 
consistent with the amplitude of temperature fluctuations seen by WMAP [H]), its actual 
value is not necessarily set by the amplitude of the fluctuations. Inflation can take place for 
field values where the shape of the potential is unrelated to the shape at its minimum. In 
addition, we make the following non-trivial assumptions about the potential near = 0: 
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Figure 1: The above figure shows an inflaton potential which can support oscillons. We are 
only interested in the shape near the bottom of the potential well. Near the minimum, the 
potential (thick black line) has to be shallower than quadratic (dotted line). The inflaton 
potential away from the minimum (dashed line) is not crucial to our analysis apart from its 
possible influence on the initial conditions. 



• V'{(f) — iv?ip < is required for some range of if for the existence of oscillons. This 
implies that A > (ignoring terms beyond ip^). Heuristically, the potential has to be 
shallower than quadratic near the minimum. 

• We assume that {\/ gY ^ 1- This is not strictly required, however it makes a semi- 
analytic analysis possible and allows for the existence of large {R ^ m~^), massive 
{Mosc ^ 'm), robust fiat-topped energy density configurations. In [IS] we provided a 
detailed analysis of the oscillon solutions and their stability in the above class of models. 
The downside of this assumption is that it makes our model somewhat special. This 
ratio appears often throughout the paper. The reader can assume a value {X/g)^ ~ 10~^ 
while reading most of the text. 

• We assume that /3 = \/X{X/g){mpi/m) ^ 1. Heuristically, /3 ~ fi/H characterizes the 
growth rate of fluctuations fi compared to the Hubble rate H. This condition on /3 
is necessary for the parametric amplification of the initial fluctuations that ultimately 
form large oscillons. We will discuss this requirement after an analysis of the linearized 
solutions. While reading through the paper, the reader can assume (3 ~ 10^. We 
keep (m/mpi) fixed in this paper, hence (3 ~ 10^ and (X/g)^ ~ 10~^ is equivalent to 
A ~ 2.5 X 10~^. When we vary parameters, we find it convenient to treat (X/g) and A 
as independent parameters instead of A and g. 
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• The coupling of the inflaton to the other fields (including the Standard Model ones) 
is assumed to be small compared to its self-interactions. This does not mean that the 
inflaton does not decay, but that significant decay happens long after the fragmentation 
of the infiaton. Since we only require infiaton decay and thermalization by the time 
BBN begins, this is easily achievable. Without this assumption, or specifying the 
nature of coupling to other fields, we cannot make any concrete statements about 
what fraction of the infiaton energy density ends up in oscillons. 

We assume that ^(0) = 0. For concreteness, from now on we will discuss the emergence 
of oscillons in the context of ([2]), with terms beyond Lp^ being explicitly set to zero. However, 
we stress that our techniques are general and can easily be applied to a wider class of models. 



3 Initial conditions and linear evolution 

In this section we follow the evolution of linear fiuctuations as they get amplified via para- 
metric resonance. For a linear analysis of the fiuctuations, it is easy to obtain the solutions 
numerically. However, we choose to provide a discussion based on approximate, analytic 
solutions to reveal the effects of different parameters on our results. 

The equation of motion for the infiaton is 

= V'{^). (3) 

We write ip{t, x) as a sum of the homogeneous piece ip{t) and fiuctuations 5({){t, x). We assume 
that {6ip{t,x.)) = where (. . .) denotes spatial averaging. Similarly the metric is FRW with 
small deviations : ^^^dx'^dx'' = {Gj^^ + 60'"') dx^'dx" where Gj^^dzf^dx" = -dt^ + a^{t)d-^- d^K 
and 5Q^y{t,yi) -C Q^y{t)- For the following, t is cosmic time in an FRW universe and x is a 
co-moving, cartesian co-ordinate on a fixed time slice. 



3.1 Homogenous background evolution 



The equations of motion for the homogeneous field are 



(4) 



where H = a{t)/a{t). We have assumed that the backreaction of the fiuctuations on the 
homogeneous equations of motion is small. That is, V" {^pji^dip"^) <^ V'{i^). After the end 
of infiation, for the potential in ([2]), the field oscillates about the minimum with a decaying 
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amplitude. At the bottom of the potential, the solution is well approximated by 



H 



where frequency of oscillation is given by cu^ ~ (m^ — ^<^^ + ...). This solution assumes that 
(f is small enough for the non-linear terms in the potential to be sub-dominant. We have 
defined a(ti) = ai = = and H{ti) = Hi. It is convenient to chose (^j and Hi (and 

hence ti) based on the structure of the instability bands for the fluctuations. The time ti is 
chosen so that before tj, there is no significant field induced amplification of fluctuations on 
scales of interest for oscillon formation. We will discuss this in a later part of this section 



[see equation (11 )]. For the moment, we assume that Hi <^m and m <^ (fi <^ nipi. 



3.2 Linear evolution of fluctuations 

The linearized equation of motion for the field fluctuations are: 

d^6ip + 3Hdt5^ + + r'iif)^ 6ip = F [Sg,,] , (6) 

where r.h.s arises due to fluctuations of the metric. Note that the linearized expression on 
the l.h.s is justified as long as V"{ip) ^ V""'((^)(5(/]^ In Fourier space we have 

d^6^k + SHdt6ipk + + V'iip)^ 5ipk = Fk [5G^,] . (7) 

The solution to the above equation depends strongly on the relationship between A;, if and 
V"{ip). When the homogeneous field is in the oscillatory regime, V^"(<^) — and F [5Q^j_v] 
are oscillatory and could give rise to parametric amplification of the fluctuations. Growth 
of 5ip was also investigated in [73] for a similar system in the late universe (without the (p^ 
term). 

Let us first consider the case where V"{ip) — m? <^ F[5Qp^v]- As a concrete example, 
consider V{ip) = m?ip^ /2. This case has been recently investigated in detail in [76] and [77]. 
On the other hand if V"{ip) — m? 2> F[5Q^j_y\ we can ignore the gravitational effects and 
concentrate on the nature of the field induced resonance H We wish to concentrate on this 
regime. 



^We thank Raphael Flauger for pointing this out 

^This is a scale and time dependent statement, and one has to look at the structure of the resonance 
bands to confirm this. We expect the metric fluctuations to play a role in setting up the initial conditions 
and also after the field resonance stops being efficient on scales k/a< yJ^Hm [751 [77] , 
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3.2.1 Initial Conditions for fluctuations 

We are interested in sub-horizon scales k/a ^ H. The initial conditions at t = depend 
on the potential V{(p) beyond <^(ti) and long wavelength {k <^ m) gravitational effects. For 
example, the amplitude of modes that never left the horizon compared to ones that did so 
during inflation can be significantly different at tj. Even if we assume that field induced 
resonance only becomes important after [this is true for our potential in ([2])], without 
specifying the detailed shape of the inflaton potential beyond (^i, it is difficult to estimate 
these effects. However, all of these effects are likely to increase the level of fluctuations 
compared to the zero point fluctuations at any given sub-horizon scale. We take a somewhat 
conservative approach in choosing our initial conditions, taking them to be consistent with 
the zero point flucutations of the field in Minkowski space. Treating these fluctuations as a 
classical, Gaussian, random field, the typical value of fluctuation on a scale k~^ is given by 
(ignoring interaction terms): 

U3/2 

eiH^, - ^ (8) 

where ul = k'^ + w? + 0[H^] (see for example [ZH]). Since we assume that Hi <^ m, we 
will ignore the Hi piece. As we will see, parametric resonance is strong enough to generate 
oscillons from these conservative initial conditions. We are treating these initial sub-horizon 
fluctuations classically, since we expect their occupation numbers to be large after they hit 
the instability band See appendix A for further details on the initial conditions and 

some justifications. 



3.2.2 Resonance in Minkowski space 

The oscillatory terms of V^"(<^) — m? in equation ([T]) can lead to parametric amplification 
of fluctuations for some selected bands of wavenumbers and values of the background field 
(f. It is these resonantly amplified modes which form oscillons. In an instability band, if we 
were to ignore expansion {a = 1,H = 0) we can treat the evolution of the fluctuations via 
standard Floquet analysis (for example, see [79]). The most unstable modes have a solution 
of the form: 

5<^fc(t)(xe'^'=*P(t), (9) 

where /i^ is the Floquet exponent and P is a periodic function. The Floquet exponent 
depends on the amplitude of the "pump" field (p as well as the wavenumber k. In figure 
^a), we show the Floquet exponent as a function of the amplitude (p and wavenumber. The 
colored regions represent the first instability band, with the color representing the real part 
of the Floquet exponent (red corresponds to a large Floquet exponent). Note that this is 
different from the usual Mathieu instability chart, where the Floquet exponents are plotted 
in terms of the resonance parameter g oc and Ak = {k/niY + 2g. Higher order, narrow 
bands [A/c < {\/gY)] exist at k > m, for example at A; ~ \/3m, beyond the right edge of 
the plot. 
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Figure 2: In figure (a), the resonance band for tlie potential V{(p) = m'^ip'^ /2 — Xif^ /4+g'^(f^ /6 
is sliown. Tlie sliaded region lias a positive Floquet exponent. Tlie tliin black-lines sliow tfie 
"patli" following by modes as the universe expands. The amplification of fluctuations can be 
estimated by integrating the Floquet exponent along these "paths". Our initial conditions 
are chosen at the dashed line. In figure (b), the amplification of a mode with {g/X)k = 0.3m 
is shown as a function of the scalefactor a. The dashed line is our analytic estimate. The 
orange lines are the two numerically integrated solutions with orthogonal initial conditions: 
6ipk = 0, 1 and dtSipk = 1,0. Also see discussion following equation (15). The waveneumber 
and field are expressed in units of the mass m. 



Under the assumption that {X/gY ^ 1 we can write down an approximate form of the 
Floquet exponent shown in figure [2|a) as 



3A 

2 \m) 



where 



- \ 2 



m I 



(10) 



3A 

5? 



is the amplitude where this band shuts off. The derivation can be found in appendix B . For 
the parameters of interest, we have checked that this expression agrees with the numerically 
calculated Floquet exponent to < 1 percent, except in a very tiny sliver at the edges of the 
Floquet band. 



8 



3.2.3 Resonance in an expanding background 



To understand parametric resonance in an expanding background, we make tlie following 
identifications: 

/c — )■ fcr, = ka~^(t), 

- - -3/2,,, (12) 



Here kp is the physical wavenumber and (fi was defined in equation (11). For future conve- 
nience, note that at = Lpi, Hi ^ ^^/X/10g^ {m/mp\) m. 



The identification ( [12| ) defines a trajectory in the kp — (p plane [thin black lines in figure 
[2|^a)]. The expression for fik{a) along such a trajectory is 




(13) 



A mode passing through the instability band can get amplified. Such modes are bounded 
by: 

< A; < 0.65 ( m. (14) 



The lower bound comes from insisting that the modes must be sub-horizon (which is a time 
dependent statement), whereas the upper bound comes from the requirement that they pass 
through the instability band [see figure |2]( a)]. 

The approximate amount of amplification undergone by a given mode is obtained by 
integrating /i^ along the corresponding trajectory in the kp — (p plane. Hence an approximate 
expression for the evolution of the amplified modes modes is 



a3/2(t) 



exp 



d(t,k) 



exp 



>, -Ma) 
dlna——— 

lJd{a,k) 



(15) 



We use the scalefactor as a time co-ordinate. The boundary of the integral d{k,a) is 
obtained from fi{k, a) = 0. We are effectively assuming that the time scale of oscillation 
(~ m~^) is much shorter than H~^. This expression should be used with caution. Its 
accuracy depends on the k mode under consideration and should always be checked with a 
direct numerical integration (also see [20]). In figure [2]^b), we compare the above expression 
with the numerical results for a particular k mode: k = 0.3{X/ g)m, with {\/gY = 0.2 and 
A = 3.13 X 10-^ H 



Ideally we would compare the numerically calculated eigenvector of the time dependent Floquet matrix 



with our approximate expression ( 15 1. However, this would require diagonalizing the Floquet matrix at every 



time step. Here the comparison is made between our expression ( 15 1 and the fundamental solutions. For the 



concerned reader, we note that our main result, regarding the number density, will depend on logk^^'^Scpk 
and not on k'^^'^Sipk. 
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3.2.4 Condition of significant amplification 



From the exponent in equation (15) it is clear that we need /iA;(a) ^ if for a significant 



amplification of the fluctuations. This means that: 



/ifc(Qj 
H 



m 



9 A;2 

10 a2 



1 - 




> 1. 



(16) 



where k = {g/Xm)k. This motivated our definition /3 = X^^"^ / g{mp\/m) in section 2j For 
appropriate values of X,g and m, we can get /3 and hence iik/H ^ 1. However, note that 
for any /3 and k, the amplification will cease eventually as a 3> 1. The amplitude can also 
stop increasing if the mode leaves the instability band, or when non-linearities cap off the 
growth. In addition, the gravitational terms from the r.h.s of equation ^ might become 
important. 

Putting all these results together, the evolution of such fluctuations after they hit the 
instability band is 



1 1 



exp 



/3fCk,a) 



(17) 



where 



f{k,a) 



2 hi 



din a 



10a2 



a3 



A;2 
a 



This integral can easily be done numerically and for a ~ a few, the integral peaks at A; ~ 0.4. 
The typical modes that are amplified have a co-moving wavenumber k < m{X/g), with a 
corresponding Floquet exponent /i ~ Q.lm{X/gY. Recall that we assumed {X/ gY ~ 0.1 and 
we require [5 = \/X(X/g){mpi/m) > 10^ for significant amplification. The Hubble parameter 
at the time of amplification is if ~ m{X/g)'^/P ~ 10~^m. Although these simple estimates 



serve as a useful guide, we use the solution presented in equation (17) for the next section 



Before we end this section we would like to visit the assumptions we made during our 
analysis. First, we have concentrated on the first instability band and have ignored the higher 
order bands (for example at kp ~ -y/Sm). The higher order bands are very narrow {Akp < 
{X/gY) and, due to expansion, modes rapidly redshift through them. The structure of these 
higher order bands is rather non-trivial. Their width in the kp—Lp plane remains small, though 
it does increase with the amplitude before shutting off. Oscillons have sizes significantly 
larger than mT^. As a result, we do not expect these high kp modes to significantly affect 
oscillon formation. Nevertheless, if modes pass through a large number of such bands or 
spend a long time in one of these narrow bands (that is, H is very small) large amplification 
is possible. These fluctuations could lead to oscillon formation if their wavelengths redshift 
by a sufficient amount. This would significantly complicate the analysis, and is not taken 
into account here. 
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Second, we ignored terms beyond <y9^ in the potential. This can change the structure of 
the resonance bands (especially for > ipi). However, one can always obtain the resonance 
bands numerically by including the extra terms and repeating the above analysis. 

Our linear analysis of the growth of fluctuations eventually breaks down as the amplitude 
of the fluctuations d^p ~ k^^'^Scpk, becomes comparable ip. At this point, we no longer have 
a homogeneous pump fleld to parametrically amplify the fluctuations. In addition, as 6ip 
approaches V"{ip)/V"'{ip), we have to take into account the interaction of different k modes. 
The modes start interacting rapidly with each other, ultimately forming oscillons. This 
process is difficult to follow analytically, however as we will see, the number of oscillons can 
be predicted based on the linear analysis. 



4 Estimating the number density of oscillons 

In this section, we estimate the number of oscillons that emerge from the breakup of the 
inflaton at the end of inflation. The idea is to use our understanding of the linear evolution 
to predict the number density of oscillons. 

Oscillons emerge from the parametrically amplifled fluctuations. As parametric resonance 
ceases to be efficient, there is a characteristic length scale where the fluctuations in the fleld 
are highly nonlinear. We claim that the number density of oscillons can be estimated using: 



2n 



(19) 



where k^i label the modes that become non-linear flrst. More explicitly, these are the modes 
for which the condition, 

k^^^6(pk ~ <f (20) 

is satisfled the earliest. We admit, that this is somewhat ad-hoc. One could have also used 
6ip ~ V"{ip)/V"'{ip) (mode-coupling) or 6(p ~ ^yV'{(p)/V"'{if') (backreaction on homoge- 



neous evolution) as a condition for obtaining a slightly different fcni. We use equation (20) 
because for the model under consideration, 5ip ^ ip happens before 5ip ~ V"{ip)lV"'{ip) 
or 5(p ~ ^JV'{ip)/V"'{ip) are satisfled. Nevertheless, we have checked that using the other 
conditions yields similar results for our estimate of the number densitjQ Heuristically, we 
are merely counting the number of large peaks in the energy density. The linear analysis 
tells us how many such peaks we should expect. 

Using k^l'^d'^k ~ (2a;fc)-i/2(A;/a)3/2exp[/3/(fc,a)] and ^ = (^iO'^/^ in (|20|) we get: 

p/2 



exp 



/3/(A^,ani) =\l^, (21) 



** A more accurate analysis would require answering a rather difficult question: Why are oscillons attractors 
in the space of many possible field configurations? 
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where ujj, = 1 + {\/ g^k?, k = m ^(\/g)k and /3 = \f\{\ / g){mpi/ m) . Different modes will 

become nonlinear at different times or scalefactors, a^\{k). We can obtain the modes that 
become non-linear first by solving for k in 



d-^aniik^i) = 0. 



(22) 



A plot of ttni VS. k is shown in figure [s^a) for different values of (3 with {X/gY = 0.2 and 
m/mpi = 5 X 10~^. The dashed line corresponds to fcni, these are modes that become non- 
linear first. Since here we have fixed m/rripi and X/g, we can treat /3 or A as an independent 
variable. They are related by A = 1.25 x lO"^*^/?^. Note that we could have equally considered 
m/rripi as the independent variable with fixed A and g. We prefer /3 because it controls how 
quickly modes get amplified. 

From figure |3]^a) we see that for a fixed (X/g), as (3 increases, k^i becomes smaller. This 
implies that the comoving number density of oscillons decreases with increasing /3. Although 
equations (21) and (22) are difficult to solve analytically, one can easily estimate k^i using 

10^ (or A = 2.5 X 10"^) we get k^i ^ OAm{X/g). This yields 



figure [3[a). For example for /3 = 
a co-moving number density of 



0.4 A 



-m 



(23) 



Note that as /3 gets large, the vs. k curve becomes exceedingly fiat near its minimum. 
As a result, k^i is not sharply defined. In this regime, our ansatz is like to fail. 

Let us understand the effects of different parameters on the estimate for the number 



density. To a good approximation, /cni obtained by solving equations (21) and (22) is given 
by: 

Ki ~ p-'/'iX/g)m. (24) 



Thus, the number density of oscillons is given by 



(25) 



where /3 = \/~X{X/ g){mpi/m). This is the main result of this section. We have checked that 
this is consistent (~ 15%) with a numerical solution of equations (21) and (22) in the range 
(A/(7)^^[5 — 10] and A[6.25 x 10~^ — 6.25 x 10~^]. We also varied m/mp\ within an order of 
magnitude of 5 x 10~® and found similar results,^ 

One might wonder if it was possible to read off k-^i from figure [2|^a) directly. The scaling 
with {X/g) comes from fc, however the dependence on /3 is somewhat difficult to see. Recall 
that (5 characterized the growth rate of fluctuations /i compared to the Hubble rate H. As 



^^Better fits can be found, but tliis is good enough for our purposes. 
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Figure 3: [3 + 1-dimensions] In (a), the solid curves show the scalefactor at which a mode 
with given k becomes non-hnear. Different curves correspond to different values of /3 = 
\/X{X/g){'mpi/m) ^ jj/H where /i is the Floquet exponent and H is the Hubble parameter. 
The dashed curve in (a) represents the modes that are the first to become non-linear. For 
the plot we chose m/m^i = 5 x 10~^ and {X/g^ = 0.2 which implies A = 1.25 x 10~^°/3^. 
Based on figure (a), in figure (b) we plot the co- moving number density of large oscillons 
as a function of /3 (or A). As expansion gets slower, the number density of large oscillons 
should decrease. Spacetime variables are expressed in units of inverse mass m~^. 



P gets larger (equivalently H gets smaller with (X/g) fixed), the low momentum modes get 
amplified before the relatively higher k modes become resonant (see shape of instability band 
in figure |2]). As a result, the non-linearity condition (20) is satisfied by longer wavelength 
modes first. 



Although plausible, equation (25) is an approximation and should be checked with de- 
tailed numerical simulations. We have not discussed some important aspects related to the 
emergence of oscillons. Non-linear interactions between oscillons is not fully understood. 
Their evolution following their emergence, as they merge and scatter is difficult to tract an- 
alytically. In particular, for a fixed (X/g), as the Hubble expansion gets smaller, (/3 3> 100, 
equivavently A ^ 10~^ or we increase m^i/m significantly), oscillon-oscillon interaction as 
well as interactions between oscillons and large non-linearities can alter the number density. 



5 Emergence of oscillons: 1 + 1-dimensions 

In the previous sections we provided an ansatz for the number density of oscillons produced 
at the end of inflation in 3 + 1-dimensions. Ideally, we would like to test this ansatz through 
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numerical evolution of the fully non-linear system on a lattice and test explore a wide range 
of parameter space. However, a large dynamic range of scales is necessary to resolve a 
significant fraction of the Hubble horizon as well as the internal structure of oscillons in an 
expanding universe. These simulations are time and memory intensive, and will be presented 
in an upcoming publication [69]. 

In this paper, we test our analysis using 1 + 1-dimensional numerical simulations. This 
significantly simplifies the numerics as well as the analytical expressions. We first reduce the 
3 + 1-dimensional analysis to 1 + 1-dimensions and then proceed towards the comparison 
with numerical simulations. 



5.1 Linear evolution and initial conditions: 1 + 1-dimensions 



The potential for the inflaton field is taken to be 



m 



1 2 A 4 02 

-IT) H w"- 



(26) 



Here h = c = 1. Note that in comparison with equation (|2]), we have scaled out m because in 
1 + 1-dimensions the field (p is dimensionless. The equations of motion in a 1 -|- 1-dimensional, 
homogeneous expanding universe are 



df^ + Hdt^ - a-'di^ + V'{^) = 



(27) 



where H = a{t)/a{t). In 1 + 1-dimensions, in the oscillatory phase, the homogeneous field 
evolves as <^(a) = ipia~^^'^ and we assume that H ^ Hia'^/"^. The choice of H amounts 
to a prescription for the expansion history since in 1 + 1-dimensions the Einstein tensor is 
identically zero. The structure of the Floquet instability band in terms of (p and kp does not 
change. Hence we take (pi = ^^/3X/5(f in analogy with the 3 + 1 dimensional case. Although 
the structure in <^ — A;p plane is unchanged, the "path" traced by a fluctuation with a given 
wavenumber does change due to the different scaling of (f and H with a. 



For a > 1, the perturbations in the field evolve as [compare with equation (17) 



k'/'SMa)-J^4^ 



1/2 



exp 



/3- 



(28) 



where = \/ 1 + C^^^, k 



m ^{g/X)k and = 1 + (10/9)A;^. We have assumed that 

Hi = ^J~Xjl^(p'{m/mp\}m = m^/l/T6{X/gy / (3. We have chosen initial conditions at a = 1 
to be identical to the 3 + 1 dimensional case. In contrast to the 3 + 1 dimensional case, we 
were able to integrate / fikit)dt analytically. 
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Figure 4: [1 + 1-dimensions] In (a), the solid curves show the scalefactor at which a mode 
with given k becomes non-hnear. Different curves correspond to different values of /3 = 
\/X{X/g){mpi/m) ^ fi/H where /i is the Floquet exponent and H is the Hubble parameter. 
The dashed curve in (a) represents the modes that are the first to become non-linear. For 
the plot we chose m/rripi = 5 x 10~^ and {X/gY = 0.2 which implies A ~ 10~^^f3'^. In figure 
(b) we plot our estimate for the co-moving number density of large oscillons as a function of 
(3 (dashed line). The orange [{X/gY = 0.1] and black [{X/gY = 0.2] data points and 2 sigma 
error bars are from our 1 + 1 dimensional numerical simulations, each point based on 10 
independent realizations of the initial conditions. The dependence on X/g is also captured 
by our estimate. The number density decreases with H. Note that A(m/mpi) is being varied 
over 2(1) orders of magnitude. Spacetime variables are expressed in unit of the inverse mass 



5.2 Number density: 1+1 dimensions 



As discussed before, oscillons form from the parametrically amplified modes. The modes 
that become non-linear earliest are the ones that satisfy the condition 6ip > first. The 
scalefactor Oni at which a given mode becomes non-linear and the co-moving wavenumber of 
the mode that becomes non-linear, kni, can be obtained from 



1 - 



6"fc A 
5feA g 



2/3' 



(29) 



dj,ani{kni) = 0. 



Repeating the procedure outlined in the 3 + 1 dimensional case, in figure |4| we show an\{k) and 
kni for different values of /3 ~ lO^^VX. We have fixed {X/gY = 0.2 and m/m^i = 5 x 10^^, 
however, the variation in {X/g) is essentially captured by our scaling of the horizontal axes. 
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Figure 5: The zoo of oscillon energy density profiles extracted from our simulations. At low 
amplitudes we see the oscillons field configuration has a "sech" profile. The larger amplitude 
objects have fiat tops (flat-top I). The widest objects (flat-top II) in the simulation show a 
slow time scale {t ^ m^^) breathing mode. Typically their amplitudes are slightly above 
the critical amplitude discussed in the text. The first two types of objects are very well fit 
by our analytic expressions for their profiles. However, we can not analytically capture the 
long term breathing mode for the widest objects. 



Unlike the 3 + 1 dimensional case, aai{k) can be written down analytically. An approximate 
expression for the co-moving number density is then given by 



2tt 



-2/5 



A m \ 



(30) 



Apart from the scaling of k by {\/g)m, the strongest dependence is due to the (3 appearing 
in the exponent of equation (28). For m/mpi[5 x 10~^ — 5 x 10~^], {X/g)~'^[5 — 10] and 
/3[30 — 300] (equivalently A = [6.25 x 10^^ — 6.25 x 10^^]), the above expression is consistent 
with the solutions of equation (29) to within ~ 15%. It is of course possible to find a better 
fit, however the above result serves as a useful guide for comparing with simulations. We 
will compare nose obtained above with the full numerical simulations. 



5.3 Numerical simulations in 1+1 dimensions 

In this section we follow the full non-linear evolution of the field numerically. Details of the 
numerical set-up and initial conditions are deferred to appendix A. First, we will review the 
properties of individual oscillons we expect to see emerging from our simulations. Second, 
we will present the numerical evolution of the field and energy density using a fiducial set 
of parameters and a particular realization of initial fluctuations. Third, we will average over 
the different realizations of the initial conditions and see how various observables evolve as a 
function of time and provide some statistical information about the individual characteristics 
of emergent oscillons. Finally, we will vary the parameters and see how the number density 
and fraction of energy density in oscillons depends on the parameters. We will see that the 
results are in good agreement (well within a factor of 2) with our analytic estimates. 
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5.3.1 Individual oscillons 



For numerical purposes, an oscillon is defined as a persistent, localized fluctuation with a local 
maximum, whose energy density at this maximum is at least 5 times the mean density. The 
width of an individual oscillon is deflned as the size of the region where the energy density 
is greater that 1/e of its value at the center of the oscillon. The energy of an individual 
oscillon is deflned as the energy enclosed within the above deflned width of the oscillon. In, 
[15] we provided an analytic solution for the oscillon proflle in an expanding universe. For 
{X/gY <^ 1, and H m, the oscillons are described by 



1 + "U 

^^''^) = ^<^^/l+.cosh[(2aA/,)x] ^^'^ 
where we have assumed X'^/g'^ ^ 1 and have ignored terms that are higher order in X/g. 



u = a/I - {a/acy 



^0 



'\{l-u), (32) 



lOg' 

co^ =m^[l- iX/gfa^] 

This is a one parameter family of solutions (once A,^' and m are specifled), whose shape 
depends on < a < a^. = a/27/160. The width is a non-monotonic function of a. It diverges 
at a — > and a — > ac- The amplitude y^o ~^ as a — )■ 0, but approaches a flnite value 
(fo — > a/QA/IO^ as a — )■ etc- The relationship between the width and the height is shown 
in the left panel of flgure |8} As a <^ 1 (same as small amplitude), we get oscillons whose 
fleld proflle is given by a sech function. However as a — )■ ac, we get a flat-topped proflle. 
Oscillons extracted from our simulations contain both forms of the solution (see flgure [S]). 
Note that a "Gaussian" proflle often used in the literature is not a good flt for the flat-topped 
oscillons. As discussed in [SnillS], our solution changes character at x* ~ a[{g/ X)H]^^ and 
becomes oscillatory in space due to expansion effects. The above solutions remains a good 
approximation as long as the width Xe < x*, else it gets stretched out by the expansion. 
For flxed set of parameters X,g and m, oscillons can exist for arbitrarily large widths (with 
amplitudes between and ■\/9A/10(7^), but expansion effects limit their sizes. The expansion 
also causes a slow loss of energy from the oscillons in the form of outgoing radiation 



We digress briefly to consider some additional properties of extremely large width oscil- 
lons, based on the linear stability analysis of |15] and [16]. In the a ^ 1 limit ("sech" like 
proflle), 3 + 1 dimensional oscillons suffer from a collapse instability when the wavelength 
of the perturbations are comparable to the width of the oscillons [13]. As a increases, mov- 
ing towards the flat-topped conflgurations, this instability disappears. In 1 + 1-dimensions, 
the collapse instability is not present. In the a ^ ac limit, extremely wide oscillons can 
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Figure 6: The above figure shows the fragmentation of the infiaton during the oscillatory 
phase of the infiaton. The large spikes in the field are oscillons. Note there characteristic 
scale fcni at which we first develop the largest fluctuations. The number density of "first 
generation" oscillons is determined by this scale. Note that the oscillons appear to get 
thinner because their physical size is fixed as the universe expands. Here the initial Hubble 
parameter k, 10~^m whereas {X/gY = 0-2 and m/rripi = 5 x 10~^. The co-moving size 
of the simulation volume is L ^ H^^ and we allow the universe to expand by a/ = 40. An 
animation of the process can be found online, 
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Figure 7: The above figure shows the fragmentation of the infiaton during the oscillatory 
phase of the infiaton. The large spikes in the energy density are oscillons. Note there 
characteristic scale fcni at which we first develop the largest inhomogeneities. The number 
density of "first generation" oscillons is determined by this scale. Note that the oscillons 
appear to get thinner because their physical size is fixed as the universe expands. Here the 
initial Hubble parameter Hi ^ 10~^m whereas {X/gY = 0-2 and m/m^\ = 5 x 10~^. The 
co-moving size of the simulation volume is L and we allow the universe to expand by 

a/ = 40. The energy density is expressed in units of m^. An animation of the process can 
be found lonlinel 
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efficiently transfer energy to {k ~ \/3m) perturbations [IH]. As a result, we expect an up- 
per limit on the width of oscillons seen in our simulations. We also note that we see some 
breathing-mode configurations of localized energy densities when a > ac, for whom we do 
not have an analytic description (see fiat-top II in figure [s]). We conjecture that they are 
likely to be bound states of two oscillons. 



5.3.2 Field and energy density evolution: single realization 

Let us now follow the evolution of the field and energy densities for a single run. In figures |6] 
and [7] we show the evolution of the field and the energy density. The values of the parameters 
were chosen to be m/mpi = 5 x 10~^, {X/gY = 0-2 and A ~ 3 x 10""^ (/3 = 50). The large 
spikes in the field and energy density are oscillons. Note that the oscillons appear to get 
thinner because their physical size is fixed as the universe expands and the horizontal axes 
are labelled in co-moving co-ordinates. The initial box size is L ^ Hr^ = 800m-^ At the 
end Lf ^ GAHj^. In the simulations H ^ Hi/y/a. 

For a ^ ^ 2.7 we do not see any significant amplification. At a ~ a^i significant 
deviations from the homogeneous energy density appear with a characteristic wavenumber 
kni- This leads to the formation of the first generation of large oscillons. Their number 



density is predicted by the formula in equation (30). A second burst is seen at 02, with 
energies smaller than the first generation ones. There could be subsequent bursts, however 
the energies of the oscillons produced tend to be significantly lower than the first burst. 
In addition, these low energy oscillons (with small amplitudes) have large widths and get 
stretched out by the Hubble expansion. 



5.3.3 Field and energy density evolution: statistics 

We ran 10 simulations with the same parameters, with different initial conditions for the 
fiuctuations drawn from a Gaussian distribution (also see appendix A). In figure [s] (right 
panel) we show a histogram of oscillon energies at a/ = 40. The distribution of energy 
densities is bimodal at late times and corresponds to the two generations of oscillons. Based 
on this histogram we call oscillons in the right lobe of the bimodal distribution, the "first" 
generation oscillons. One might wonder if the lumps we are seeing are indeed oscillons. We 
can get rid of of the false positives by comparing the height and width of the energy density 
of the lumps with the theoretical relationship (see the left panel in figure [s]). As seen in 
figure [7| at late times, all large lumps that remain are oscillons. Again, note the bimodal 
distribution of widths and heights of the oscillons corresponding to the two generations. 



In figure [9] (right panel) we show the number density of oscillons (total and first genera- 
tion) as a function of time. Large oscillons are produced at ~ a^i and thereafter freeze out 



with the expansion. The dashed line is our analytic estimate in (30). The results of fraction 



of energy density in oscillons (after averaging the runs over 10 realizations of the initial 



20 




Figure 8: [1 + 1-dimensions] In the left panel, we show the non-monotonic relationship 
between the width and core amplitude of the oscillon energy densities based on our analytic 
solution (solid black line). The orange points are from the energy density configurations 
flagged as oscillons in our simulations. At large core amplitude, we approach a fiat-topped 
profile. In the right panel, a histogram of the energies of oscillons is shown [A ~ 3 x 
10~^, {X/gY = 0.2,m/mpi = 5 x 10~^ with 10 different initial condition realizations]. Two 
distinct oscillon populations are clearly visible, with the larger energy ones appear first in 
the simulations. All variables are expressed in units of appropriate powers of mass m. An 
animation of the process can be found jonline 
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Figure 9: [1 + 1-dimensions] The left panel shows the number density of oscillons (total and 
first generation) as a function of the scalefactor. Note that they are produced at a^i ~ 2.7 and 
thereafter freeze out. The dashed line is our analytic estimate in (30). The resulting energy 
density in oscillons after averaging the runs over 10 realizations of the initial conditions is 
shown on the left. The 2 sigma error bars quantify the variations between different runs. 
Note that most of the energy density is in the first generation of oscillons. The values of 
the parameters used were A ~ 3 x 10^'^ (or /3 = 50), {X/gY = 0-2 and m/mpi = 5 x 10^^. 
Number density is expressed in units of the mass m. 



conditions) is shown in figure |9] (left panel). The 2 sigma error bars quantify the variations 
between different runs. Note that most of the energy density is in the first generation of 
oscillons. The values of the parameters used were A ~ 3 x 10~^ (or /? = 50), {\/g)^ = 0.2 
and m/mpi = 5 x 10^®. 



5.3.4 Field and energy density evolution: parameter dependence 

We will now vary the parameters of the model and see if our results match the semianalytic 
estimate. Figure |4] (b) shows the final number density of oscillons at a/ = 40 as a function 
of (3 for (X/g)^ = O.l(black) and {X/gY = 0.2 (orange). The dashed curve is our analytic 
estimate. It is somewhat remarkable, that in spite of the extremely non-linear dynamics 
of oscillon formation, our analysis gets the number density to well within a factor of 2 and 
correctly captures the variation with parameters. 

At small P < 20, parametric resonance stops being effective and we get fewer and smaller 
(in amplitude) oscillons. At large (3 we notice a systematic deviation from our analytic 
estimate. As seen in figure |4](b), the analytic result over-estimates the number of oscillons 
at large /3. However, at large /3, the rate of expansion is slow. As a result, in our simulations 
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we get a large number of oscillon-oscillon interactions, which leads to some of them being 
disrupted. This was not taken into account in our estimate. Oscillons can merge, scatter 
off each other or be disrupted, depending on their amplitudes and phases. Our exploratory 
investigation with oscillon collisions indicates that the outcome of the collision depends 
on the phase, amplitude as well as shape of the oscillons (also see [38] )■ The collisions 
range between almost completely elastic collisions to highly inelastic ones (though we are 
almost always left with one or more oscillons). The details of oscillon-oscillon interaction 
or interactions with other large non-linearities is not completely understood and certainly 
warrants further investigation. 

A second effect is that as fcni gets smaller, we get very wide oscillons. These objects suffer 
from a Floquet instability at ~ \/3m HE] and can potentially destabilize the oscillons. 
We also note that in our simulations, the largest width objects (see flat-top II in figure 
[5]) have an energy density that stays localized but exhibits a long time scale {t ^ m^^), 
breathing mode which do not fully understand. 



The oscillons tend to dominate the energy density of the field. In figure 10 (right), we 
quantify this statement. As expected, at very small /3 we cannot generate oscillons efficiently, 
so this fraction must increase as a function of p. For the parameter range considered, oscillons 
take up more that 75 — 90% of the energy density of the field. Note that for numerical 
purposes we defined the oscillon energy as the energy contained within a width where the 
energy density is above 1/e of its core value. Thus, we we are ignoring the energy density in 
the tails. 

The energy of typical oscillons (at a fixed X/g) increases with increasing /3. That is, as the 
Hubble parameter gets smaller, we get larger, more massive oscillons without a significant 
change in the fraction of energy density in oscillons. The mean energy of first generation os- 



cillons is shown in figure 10 (left). The two colors represent different values of {X/ g)"^ [0.1, 0.2]. 
The mean and two sigma error bars are obtained from a sample containing a few thousand 
oscillons (100 realizations, 10 for each value of /3) . 



6 Discussion 



In this paper we have investigated the post inflationary, emergence of oscillons in a class of 
single fleld inflaton models. We gave analytic results for 3-1-1 and 1 -|- 1-dimensional cases 
and numerical results for the 1 -|- 1-dimensional case. Starting from zero point fluctuations of 
the inflaton during its oscillatory regime, we provided an (approximate) analytic description 
of the linear evolution of the fluctuations and provided a condition necessary for signiflcant 
ampliflcation of the fluctuations. Using this linear analysis we calculated the characteristic 
scale which is the flrst to become non-linear. We hypothesized that it is this scale that 
determines the number density of oscillons. We checked our analysis in detail with 1-1-1 
dimensional simulations in an expanding universe, varying different parameters over an order 
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Figure 10: [1+1- dimensions] The left panel shows the percentage of energy density in oscillons 
at the end of our simulations (a/ = 40) for different values of /3 = \/X{X/ g){mp\/m). The 
orange and black curves with 2 sigma error bars correspond to simulations with parameters 
{\/gY = 0-1 cind 0.2 respectively. Note that for our fiducial value of m/mp\ = 5 x 10~^, we 
have A ~ 10~^°/3^. The fraction of energy density in oscillons decreases at low (3 (large H) 
because resonance is not efficient enough to significantly amplify initial fiuctuations. The 
right panel shows the mean energy of the first generation oscillons. Note that it increases 
with increasing /3 (decreasing H). This is consistent with the idea that at large /3 oscillons 
form from longer wavelength fiuctuations, thus yielding fewer, larger oscillons. Energy is 
expressed in units of the mass m. 
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of magnitude. Here we found agreement between our analytical and numerical results to well 
within a factor of ~ 2. The number density decreased and the size of oscillons increased 
with decreasing Hubble (all other parameters fixed). A detailed analysis revealed that the 
individual characteristics of the oscillons (in particular the natural emergence of flat-top os- 
cillons) extracted from our simulations were in excellent agreement with the analytic results 
of our previous paper We found that the fraction of energy density in oscillons as the 
parameters were varied was 75% — 90%. We also pointed out some interesting phenomenon 
seen in the simulations which we cannot completely account for quantitatively. These in- 
cluded the production of oscillons in more than one burst (however, see [80]), oscillon-oscillon 
interactions at slow Hubble rates and the slow breathing modes of extremely wide, fiat-top 
oscillons. 

There are many ways in which our analysis could be extended. In particular, we ignored 
the effects of higher order resonance bands. Although reasonable for the model under con- 
sideration, this need not be true for other models and should be included in the analysis. 
The effect of coupling to other fields needs to be investigated PU] . A more careful analysis of 
the initial conditions and the evolution of fluctuations including gravitational perturbations 
is also needed. Most importantly, our numerical analysis was done in 1 + 1-dimensions. 
In an upcoming publication [69], we will provide detailed numerical results for the 3 + 1 
dimensional case and compare it to the analytical estimates provided in this paper. 

In summary, to understand the cosmological consequences of oscillons, it is important to 
have a prediction for their number densities as well as their individual characteristics. For 
the model under consideration, we provided both in terms of the parameters in the infiaton 
Lagrangian. The techniques developed here should be directly applicable in a broader class 
of models. 
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Appendix A: Numerical set-up and initial conditions 



In this appendix, we outline our numerical set-up for the 1 + 1 dimensional simulations. For 
numerical purposes, we find it convenient to work with dimensionless spacetime variables 
mx^ as well as the scaled field (pp = X^^'^ip. We will work with conformal time dr] = a~^dt. 
Under these changes, the equations of motion now become: 



1 2 1 4^ 1 6 (33) 

where the 'prime' stands for derivatives with respect to conformal time and ( = X/g. To 
avoid clutter, we will drop the p subscript. Using conformal time is particularly convenient 
for numerical purposes since it gets rid of the linear derivative term. 

We will discretize the above equations in space and conformal time. We will denote 
ip{x,ri) = ip{iAx, jArj) = ipij and 0(77) = a{jAri) = aj. The symmetric space and time 
derivatives become: 



2Ar] 



(34) 

II _ V?ij+i - 2v9jj + V9i,j-i ^ ' 



Arj 

The evolution equations are given by: 

ifi^j+i = 2Lpij - (pij_i + [Pi-i^j + Vi+i^j - 2ipij] - {AT]fa'^jV'{Lpij), (35) 

where s = Arj/Ax. To evolve the system forward in time we need y^i.o and ipi^i. (pi^ is 
the initial value of the field tpi^ = ip{0,iAx) whereas ipi^i be constructed from (pi^ and 
f'io = f'iOjiAx) as follows 

= <fi,o + fifl^V + - 2v5i,o + V^j-i.o) - A7]W {(pifi) (36) 

We employ periodic boundary conditions. The initial field value and it's derivative are 
specified using (see fSU]): 

-"/^^^ (37) 

. ^/2 / ^ ^ 

^ -N/2+1 ' 
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where kn = 2nL and = y 1 + (2 sin ^"^^ / Ax) . Here, a„ are complex numbers whose 
phases drawn at random from (0, 2tt] whereas their amplitudes are drawn from the Gaussian 
distribution with variance (|anP) = A/2, consistent with the description of a massive field 
in it's ground state. A number of comments are in order. We have expressed /c„ and Un 
in units of m. We have ignored the interaction terms in specification of the ground state. 
We do not know of any prescription where these can be taken into account for the zero 
point fluctuations. The appearance of A in the amplitude of the fluctuation is somewhat 
misleading since we did not take into account the interaction terms. Here it appears simply 
because we are working with the scaled version of the field cpp = \/X(p. We chose to treat the 
evolution classically, with the anticipation that the occupation number per mode will grow 
rapidly as they undergo parametric resonance (see for example [H]). 

In addition, we have specified the initial conditions in Minkowski space. This is rea- 
sonable since, the scales of interest are much smaller than H^^. However, this significantly 
underestimates the fluctuations close to the scale of the Horizon. The nature of the spectrum 
there depends on the details of the evolution of the field before tj and the shape of the po- 
tential beyond (pi. Since we have not specified these in this paper, we take this conservative 
prescription of the initial conditions as our starting point. Another approximation is present 
at /c 3> m. As modes redshift we need to continuously re-populate the high k modes on the 
lattice. Since these high k modes do not undergo efficient parametric resonance (based on 
the structure of the Floquet diagram), we do not re-populate these modes. One has to be 
somewhat cautious here since some oscillons have discrete narrow band instabilities in the 
k ^ m region H5] . This will not be captured by the simulation. 

Now, we need to choose the spatial and temporal resolution. Oscillons maintain a fixed 
physical size as the universe expands. This means that we need to improve our spatial 
resolution as the universe expands. We do so by doubling the grid points every time the 
universe expands by a factor of 2. We make sure At] < Ax is always satisfied, by refining 
the time step along with the spatial resolution. The approach is similar to the one used in 
[50] . We interpolate the field between adjacent grid points to improve the spatial resolution, 
but use the equations of motion themselves (with assymetric time steps) for improving the 
temporal resolution. We need to make sure that we capture the fastest spatial as well as 
temporal oscillations of the system, both of which should be (at least) smaller that the 
inverse mass scale of the problem. We start with an initial grid spacing of Axq = 0.5. The 
initial time step: A?7o = 0.2Axo. By construction we cannot resolve spatial structures much 
smaller than 0.5. The size of the initial box is Lq = 800. We allow the simulation to run till 
a/ = 40. Our overall energy conservation (including the effects of expansion) is at the level 
of one part in 10^. Almost all field configurations that are flagged as oscillons are resolved 
with > 20 grid points. We have varied the spatial and temporal resolution to make sure 
that no significant qualitative differences are seen. Nevertheless, the final positions of the 
oscillons do shift as the resolution is varied. 

To evolve a{r]) we need H at each time step: a^+i = aj (1 + ajHAr]). This can be obtained 
using the Friedmann equation H = l/(3mp[)(p) where (p) stands for the spatially averaged 
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energy density. In practice, we find that H oc a~^/^ reproduces the expansion history to a 
percent level accuracy. As mentioned in the main body of the text, in 1 + 1-dimensions this 
amounts to a prescription for the evolution of the background since in 1 + 1-dimensions, the 
Einstein tensor is identically zero. 



Appendix B: Floquet exponent 



In this section we will derive the Floquet exponent in equation (10). The approach is similar 



to the one adopted in |8T] for a parametrically excited harmonic oscillatoi[^^| We start 
with the equation of motion for the homogeneous field and linear fiuctuations around the 
homogeneous solution: 



Q 

d^(f + rr?Lp - \lp^ + = 0, 



2 o \ -3 I -4 

m — ?>Xip + —7^^ 



(38) 



6ip = 0. 



It is convenient to define a dimensionless, scaled version of the field 

= yx-. 



m 



(39) 



and dimensionless spacetime variables x'^ — )• mx'^. In terms of these variables we have 



dU6 + 



-V^ + 1 - 302 + 5 



(40) 



0. 



We will be working under the assumption = {X/gY ^ 1- Under this assumption, the 
homogeneous background equation has a solution of the form: 



0(t) ^ e^o COS ujt + 0[e% 
where $o can be of order unity. The frequency of oscillation is 



w ^ 1 — e 



2/3^2 5 ^4 



The equation of motion for the perturbation 50 (in Fourier space) is 

9^50 + e -e^ (3$^ cos^ ut - 5$^ cos^ ut)] 6(j) = 0. 
WWe thank Raphael Flaugher for pointing us to this reference. 



(41) 



(42) 



(43) 
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We can rewrite the above equation as: 

dlSct) + [nl + {(3 cos 2ut + 7 cos Aut)] 5<p = 0, (44) 



where 



7 = "t^K 



2^2 

5 



We shall look for solutions of the form 

5(/)(t) = [a„,(t) cosrafifct + 6„(t) sincosnfifct] (46) 

n=l,3... 

where a„ and 6„ are slowly varying compared to the oscillatory terms. The even n terms 
are decoupled from the odd-terms and can be set to zero. Plugging this form of the solution 
into equation (44), dropping the second time derivatives (a„ and &„) and only keeping terms 
up to order e^, we get: 



^J^f (3 + 2(a;2 - nl) \ f 

h Aio[ -P + 2(^2 _ ^2) Q 11^^ 



0[e^] (47) 



The system can be easily diagonalized to obtain the following solution, conveniently 
expressed in its Eigen-basis: 

where the eigenvalue 



f^k = - Aiu^ - (49) 

is the desired Floquet exponent. Note that after an initial transient, 50 oc e'^'=* (ignoring the 
oscillatory piece). The higher harmonic terms are higher order in e. 

To lowest order in e, the Floquet exponent is given by 



2 2 "V 3 



Reverting back to the original, unsealed variables used in the main body of the text [cp 
^me$o/V^j cos cut and k — )■ k/m), we have 



k /3A /^\2 / /^x2 



where ifi = ^^/3X/5g^m. 
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